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I.  INTRODUCTION 


The  associated  Legendre  functions  on  the  cut  for  integral  order  and 
complex  degree  are  required  for  several  important  mechanical  applications. 

One  area  of  particular  interest  is  in  the  calculation  of  biharmonic  functions 
for  the  hollow  cone.  Such  functions  appropriate  for  the  cone  satisfy 


V 
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X  * 


0 


and  come  in  the  form 


X=[AP™(cos0)  +  BQ^(cosQ)  +  Cp“_2(cos0)  +  Dq“_2(cos0)]  RV  sin  m<j>. 


These  functions  are  required  for  stress  analysis  in  the  minihat  gage,  which  is 
a  strain-type  pressure  transducer,  and  in  cone-plate  rheometers. 

We  have  developed  a  subroutine  for  computing 


P^(cosG)  and  Q™(cos0) 


for 


O<0<  90°  , 

integer  m  such  that  0<m<12  , 
and  v  such  that  Re(v)>  -  V2  • 


For  moderate  values  of  |v|,  we  calculate  P^(cos0)  by  series  and  recurrence 
formulas.  On  the  other  hand,  Q™(cos0)  is  computed  through  the  use  of  a  Gauss 
continued  fraction  and  Wronskian  relations.  For  small  angles,  the  logarithmic 
solution  is  used  to  calculate  Q^(cos0).  Additional  analysis  will  be 

3Qm(cos0)  3pm(cos0) 

necessary  for  the  programming  of  - - -  and  - r -  .  Asymptotic 

<3  v  oV 

expansions  are  required  when  |v|  is  large. 

Programming  was  done  on  a  CDC  7600  computer  in  FORTRAN  V.  CDC  double 
precision  arithmetic  was  used  to  get  the  greatest  accuracy  possible.  Special 
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mul tipLe  precision  algorithms  were  avoided  so  that  the  subroutine  would  be 
more  compatible  for  use  by  an  engineering  laboratory. 


II.  ASSOCIATED  LEGENDRE  FUNCTIONS  OF  THE  FIRST  KIND 
A.  Hypergeometric  Series 


The  Legendre  functions  of  the  first  kind  can  be  accurately  calculated  for 
moderate  |v|  by1 


Pv  (cos^)  mj  (tan  0/2)  2 ^ (-v , v+1 ; 1+m; ( sin0/2) ^ ) 


(1) 


and  the  reflection  formula^ 


P  (cosO) 
v 


^  r(v+m+1)  r“m/  nv 

(  1}  T(v-m+I )  Pv  (cos0)  » 


(2) 


where 


f  (  v+m+ 1 ) 

p(v_m+1)  =  (v+l-m)(v+2-m)  . . .( v-hn) 


(2a) 


The  hypergeometric  series,  for  small  |v|,  (|v|<2),  can  be  calculated 
directly  from  the  series^ 


,F  (a,b;c;z)  =  1  +  —  fy  +  £.(.a+l>b(b+l) 

-  1  cl!  c(c+l) 


z 

l\ 


+ 


(3) 


For  larger  values  of  |v|  ,  the  series  becomes  slowly  convergent 

round-off  error.  In  this  case,  one  can  alter  the  parameters  and 
recurrence  formula  to  obtain  the  desired  results. 


with  serious 
then  use  a 


A.  Erdelyi,  et.  at .,  Higher  Transcendental  Functions,  Vol.  1 ,  Bateman 
Manuscript  Project ,  McGrow  Holt  Book  Co.,  Inc.,  'Jeo  York,  TWSZ . 

Wj.  Magnus  and  F.  Oberhettinger,  Formulas  and  Theorems  for  the  Special 
Functions  of  Mathematical  Physics,  TdKeUsed  FuFTtsFiTng  Co.,  ’New  ¥ork."~l 94 9 . 

M.  Abramoioitz  and  J.A .  Stegun,  Handbook  of  Mathematical  Functions  with 
Formulas,  Graphs,  and  Mathemat^a l^T^Tes,  Wainumdr~Buredu  of  SEdhaarTs, 
Washington,  D.C.,  No~.  65~,  AppTTed  Mathematics  Series,  June  1964. 
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B.  Gauss  Contiguous  Formula 


One  possible  method  to  improve  the  convergence  of  the  hyp ergeome trie 
series  would  be  to  increase  m  to  a  high  value  and  then  use  a  Legendre 
recurrence  formula  to  obtain  the  desired  answer.  Unfortunately,  either 
serious  underflow  or  overflow  problems  would  begin  to  occur  in  the  calculation 
of  the  initial  values  for  the  Legendre  recurrence  formula,  as  found  by  Smith, 
Olver,  and  Lozier.^  Occasionally  both  underflow  and  overflow  occur  in  the 
same  sequence  of  calculations.  The  necessary  recurrence  formula,^ 


P^+1(cos0)  +  2mcot8P^(cos8)  +  (v+m) ( v-m+1 )P™  ^(cosQ)  =  0,  (4) 

also  contains  both  linear  and  quadratic  coefficients  which  would  lead  to  rapid 
variations  of  the  Legendre  functions  in  a  recurrence  chain. 

In  order  to  avoid  the  difficulties  in  dealing  with  the  Legendre  function, 
we  deal  directly  with  the  hypergeometric  series  using  the  Gauss  contiguous 
formula.  In  the  required  Gauss  contiguous  formula,* 


c(c-l  )(z-l  )F(a,b;c-l  ;z)-fc  [c-l-(2c-a-b-l)z]  F(a  ,b;c;  z) 
+(c-a)(c-b)z  F(a,b;c+1 ; z)=0  , 


(5) 


the  coefficients  of  the  hypergeometric  functions  are  of  the  same  degree  in  the 
parameter  c.  Therefore,  the  values  of  the  hypergeometric  functions  vary 
slowly  as  the  parameters  change  in  the  recurrence  chain.  Also,  the  stability 
criterion  for  the  Gauss  contiguous  formula  can  be  verified  by  elementary 
methods.  Finally,  through  the  use  of  the  hypergeometric  functions  as  the 
basis  of  our  analysis,  we  gain  an  additional  parameter  leading  to  greater 
flexibility  in  deriving  algorithms. 

Since  the  Gauss  contiguous  formula  appears  to  be  superior  for  the 
intended  work,  it  is  used  in  the  program.  We  increase  c  to  a  sufficiently 
large  value  and  then  increment  downwards  with  the  Gauss  contiguous  formula  to 
the  desired  result.  Using  this  technique  has  produced  excellent  results, 
provided  that  |v|  is  not  too  large. 


J*M.  Smith,  E\  V.J.  Olver  and  D.W .  Lozier,  ”Legendre  Polynomials,”  ACM 
Transactions  on  Mathematical  Software,  7,  94-95,  1981 . 

A.S.  Elder,  ” Formulas  for  Calculating  Bessel  Functions  of  Integral  Order 
and  Complex  Argument,”  USA  Ballistic  Research  Laboratory ,  Report  No .  BRL 
TR-1423,  1968  AD  # 680209 . 
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III.  ASSOCIATED  LEGENDRE  FUNCTIONS  OF  THE  SECOND  KIND 


A.  Gauss  Continued  Fraction  Method 

In  attempting  to  calculate  Q™(cos0),  no  series  was  found  which  converges 
well  for  small  v  and  0.  Instead,  an  approach  previously  used  to  calculate 
Bessel  functions  of  the  second  kind  was  used.  This  approach  utilizes  the 
Gauss  continued  fraction  and  Wronskian  relations.  It  was  found  to  be  highly 
advantageous  to  leave  the  cut  for  the  initial  calculations  and  work  in  the 
complex  plane.  However,  Q™(x) ,  as  defined  on  the  cut  -l<x<+I,  is  not  analytic 
when  x=z,  a  complex  number.  We  have  the  formula* 

Q™Cx)  -  l/2  e-1"’  [e-  V2  i¥(£<*-i(»]  ,  (6) 


where  x+10,  x-iO  indicate  limits  as  the  complex  variable  z  approaches  the 
cut.  If  v  —  a+i3,  and  3  is  large, 


Q™(cos0+iO)  is  0(e 


and 


Q™(cos0-iO)  is  0(e+^Q),  where  x  =  cos0. 


This  behavior  corresponds  to  that  of  Hankel  functions  in  the  theory  of  Bessel 
functions • 


In  the  complex  plane,  we  have 


1 


e-imu  Qm  (cog  Q+iQ  >  r  (  v+3  /  2  >  =  V2  r  (  v_hn+ 1 )  eV2  mTI  (  2  s  in  0  )ra 


(7) 


—i0( 1+v+m) .  .0/0  “2i0v 

x  e  y  F  (  V; 2  >  1+v-Hn ;  v  f  3  /  2 ;  e  )  . 


Replacing  v  by  v-1  and  then  dividing  the  result  into  the  previous  equation. 


Qv(cos9+iO)  (v+u)  -i.0  F(  l/2+  m,  1+vfm; v+3/2;  e  J 

Rm(cos0*fiO)  =  -  - - e - - 

Qm  (cos0+iO)  (v+  V?  )  p(  V2  +ra>  v+m;vf  V2  ) 
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b  e 
c 


(9) 


R™(cosO+iO) 


-10 


F(a,b+1;  c+1 ;  z) 

F(a,  b;  c;  z) 


This  can  be  calculated  by  using  the  Gauss  continued  fraction 


6 


F(a,b+l;c+l;z)  =  1 _ 

F(a,b;c;z)  a(c-b)z 

1  ~  c(c-t-l) _ 

1  -  (b+l)(c-a+l)z 

(c+l)(c+2) _ 

1  -  (a+l)(c-b+l)z 

(c+2)(c+3) _ 

1  “  (b+2)(c-a+2)z 

(c+3)(c+4)  , 

1  -  ... 


The  Gauss  continued  fraction  converges  everywhere  in  the  complex 
on  the  real  axis  from  +1  to  +  00  and  at  isolated  zeros  of  F(a,b;c 

In  order  to  obtain  Qm(cos0+iO)  from  its  continued  fraction, 

V  o 

the  following  Wronskian  relations  : 


dQm(z) 


dP^z) 

-  Q  (z)  — i - 

xvv  7  dz 


and 


W{p;m(z),  Q™ ( z ) } 


iTTm 

e 


(l-z2) 


6 


(10) 


plane,  except 
z)  . 

we  utilize 


(id 


H.S.  Wall,  Anal^tVo  Theory  of_  Continued  Fnaotions,  d.  Van  Nostvand  Company 
Ino .,  Ne<J  York,  1948. 


Equating  these,  we  have 


d(f(z)  dP^(z) 

Pv  (z)  ~di - Qv(z)  ~di - 


nrm 


<l-z2>  ' 


Now  using  the  derivative  relations 


1 


.  dP^(z) 

( z;2— 1 )  — ^ -  =  v  z  P>)  -  (v-m)  P^(z) 


and 


.  2  ..  «>> 
(z  -1) 


v  z  Q™(z) 


dz 


(v+m)  Q^Cz)  , 


we  find 


(12) 


(13) 


(v-hn)  P-“(z)  (z)  -  (v-m)  (z)  q”(z)  =  e1™  .  (14) 

Dividing  both  sides  by  Q™_^(z)  and  replacing  z  by  cos0+iO  ,  we  get 

inm 

(v-hn)  P_ra(cos0+iO)  -  (v-m)  P  m1  (cos0+iO)  R™(cos0+iO)  = - - - . 

v  v  1  v  ^(cose+iO) 

(15) 


Replacing  v  by  v+1  and  then  solving  for 


Q™(cos0+iO) 


> 


Q™(cosQ+iO) 


inm 


_ e _ 

[(v+m+l)P^1(cos0+iO)-(v-m+l)P^n(cos0+iO)R™+1(cos0+iO)  ] 


(16) 


Since 


1 


10 


(17) 


Pv  (cos0+iO) 


P  1  (cos  0  ) 

V 


Q^(cos0+iO) 


V?  inur 

e  ^ 


[  ( v-hn+1  )P^^(cos0  )-(v-m-fl  )P^m(cos0  (cos0+iO)  ] 

(18) 


Now  we  can  use  the  relation 


1 


„ v  3/2  imu  ,  „  ,  . ^ N  .  ,  . _  m ,  ~  N 

Q^(cos0)  =  e  Q^(cos0+iO)  +  iir/2  P^(cos0) 


(19) 


to  finally  get 


Q™(cos0) 


- ( — l — - +  Itt/2  Pm(cos0). 

[(v+m+l)P^^  (cos0)-(v-m+l  )P^m(cos0)R™+^  (cos0+iO)  ] 


(20) 


Results  from  programming  the  above  method  indicate  that  it  works 
exceptionally  well.  Unfortunately,  this  method  will  likely  run  into 
difficulties  for  large  |v|.  Even  though  the  Gauss  continued  fraction  seems  to 
converge  more  rapidly  with  increasing  |v|,  difficulties  may  be  encountered 
because  of  losses  in  accuracy  with  P^(cos0).  For  large  |v|,  it  will  be 
necessary  to  use  asymptotic  expansions  to  calculate  Q™(cos0)  and  P™(cos0). 

Another  region  where  trouble  occurs  is  at  small  angles  (0<1°).  The 
complex  Gauss  continued  fraction  algorithm  is  still  very  accurate  at  small 
angles,  but  the  continued  fraction  converges  at  a  very  slow  rate.  Sometimes 
several  thousand  terms  are  required  for  convergence,  which  in  turn  uses  an 
inordinate  amount  of  computer  time.  This  problem  occurs  because  the  argument 
approaches  the  cut  of  the  Gauss  continued  fraction  at  +1.  To  avoid  this 
problem,  we  use  the  logarithmic  solution  to  calculate  Q^(cos0)  for  small 
angles  . 

A  similar  method  for  computing  Qm(cos0)  has  also  been  developed  using 

v  ni 

real  analysis  only.  Since  no  simple  series  formula  for  Q^(cos0)  could  be 

found,  we  concentrated  on  Pm(-cos0)  instead.  Pm(-cos0)  is  related 

v  v 

to  Q™(cos0)  and  p™(Cos0)  by1 
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IT 


q“(cos9) 


2sin[7r(v+m)] 


{cos  [ir(v+m)]P^(cos0)  -  P™(-cos0)}, 


0<  0  <TT  /  2  • 


(21) 


Similar  computational  difficulties  exist  for  the  computation  of  P^(-cos0)  as 
Q™(cos0).  Using  the  same  equations  as  in  calculating  P™(cos0),  we  can 
derive  the  following  continued  fraction 


p“(-cose) 

pv-i(-cos0) 


(v+m)  F(-v>v+l;  1+m;  (cos9/2)2) 

(V_m)  F(— v+1 ,  v;  1+m;  (cos0/2)2) 


Letting  a=-v ,  b=v , 


c=m, 


z=(cos0/ 2) 


2 


then 


(22) 


m .  .  F(a,b+l;c+l;z) 

'  F(a+1 ,b;c+l;z) 


F(atb+l;e+l;z)  F(b,a;c;z) 
F(ajb;c;z)  F(b ,a+l ;c+l ; z) 


(23) 


This  can  be  calculated  by  either  two  Gauss  continued  fractions  or  by  two 
continued  fractions  of  Frank^  which  are  of  the  form 


F(a,b;c;z) 

F(a,b+l;c+l;z) 


(Equation  is  incomplete,  continued  on  next  page) 


2  E.  Frank ,  "A  New  Class  of  Continued  Fraction  Expansions  for  the  Ratios  of 
Hyper geometric  Functions , "  Transactions  of  the  American  Mathematical 
Society ,  81_  ,  454,  1956. 
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1 


a(c-b)z 

c(c+l) _ 

(a-b)z  (a+l)(c+l-b)z 

(c+1)  1  (c+l)(c+2) _ 

(a+l-b)z  (a+2)(c+2-b)z 

(c+2)  (c+2)(c+3) 

(a+2-b)z  ,  , 

i  i  “  •  •  •  • 

(c+3)  (24) 


The  continued  fraction  of  Frank  has  its  cut  on  the  unit  circle.  It  has  proved 
to  be  of  little  value  here  because  of  its  much  slower  convergence  rate  as 
compared  to  the  Gauss  continued  fraction.  This  is  especially  true  for  small 
angles,  which  effectively  limit  the  use  of  the  Frank  continued  fraction  to 
angles  greater  than  30°. 

Utilizing  similar  Wronskian  relations  as  in  the  complex  analysis,  we 
proceeded  in  a  similar  manner,  with  the  exception  of  replacing  Q^(cos0)  by  Eq. 
(21).  The  final  result  of  the  real  analysis  is  then 


Q^(cos0)  =  P^(cos0)  cot  [tt( v*fm)  ]  + 


(~l)m _ 

(v+l+m)  [pJ^CcosQ)  +  P^CcosejS^Cz)] 


(25) 


This  result  can  easily  be  seen  to  be  inferior  to  the  complex  analysis 
since  it  is  undefined  for  integral  v  .  Also,  it  is  necessary  to  calculate  two 
continued  fractions  which  have  relatively  slow  convergence  rates  which  worsen 
with  increasing  |v|  .  Therefore,  only  the  complex  algorithm  is  used  in  the 
final  programming.  The  real  method  will  be  useful  in  helping  to  check  and 
compare  the  accuracy  of  the  complex  method  for  calculating  Q^(cos0). 

B.  The  Logarithmic  Solution 

As  mentioned  previously,  the  logarithmic  solution  is  used  for  the 
calculation  of  Q^(cos0)  at  small  angles  (0<i°)#  In  this  range,  the 
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logarithmic  solution  is  as  accurate  as  the  Gauss  continued  fraction  method, 
but  requires  much  less  computer  time  to  obtain  results.  From  references  [1] 
and  [2],  the  following  formulas  can  be  derived. 


Q”(cos0)  =1/2P™(cos0)[log(l/tan2(0/2))  -  2y  -  2*(v+l)  +  Y  - T  ] 

v  v  (v+s)(v+l-s) 


+  2^r~  I  (-v2)(l-v2)(4-v2)...((r-l)2-v2)(xH-r)-^IflUI 


r=0 


r! 


x  (~l)r(sin0/2)2r 


m  <x>  9  9  9  9 

(sin0/2)m  y  (-v  )(l-v  )...((m+£-l)  -v  )(m+v+Ji) 


2v 


x  (1  +I  +  I+  ...  +  |)(sin0/2) 


21 


(~l)m 

2v 


(  v-m+1 )  ( v-m+2)  . .  •  (  v-hn)  ( tan0/  2  Jm 


x 


00 


I 

r=0 


(~v) 2(l-v2)...((r-l) 2-v2) ( v+r) 

v!  (r+m)! 


» 


+ 


1 

m+r 
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(26) 
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14 


For  m=0,  bat  v  not  an  integer. 


Q°(cos9)  =  l/2  P°(cos0)  [log(cot2(0/2))  -  2y  -  2\Kv+l)] 


00 


(-v2)(l-v2) ...f(£-l)2-v2)(£+v) 

(£!)2 


(l+  \  \  +...|)(sin0/2)2£  . 


(26a) 


For  m=0  and  v=n  ,  an  integer, 


Q°(cos9)  =V2P°(cos0)  [log(cot2(0/2))  - 


2(1+  \  +  j  +• . .  1/n)  ] 


+ 


n  £  I  i 

v  (-l)*(n+£)!  (1+  2  +  3. 

£=1  (£!)2(n-£) ! 


.£) 


(sin0/ 2) 


2£ 


(26b) 


IV.  CONCLUSIONS 

We  have  developed  an  accurate  and  efficient  method  of  calculating  the 
associated  Legendre  functions  for  integral  order  and  complex  degree.  In 
general,  the  accuracy  of  P™(cos0)  and  Q™(cos0)  was  approximately  equal  to  the 
maximum  attainable  with  double  precision  arithmetic.  By  comparisons  with 
special  cases  (P_  1/  +  ^(cos0),  Qy(cos0)  for  real  v,  Legendre  Polynomials)  and 
alternate  means  of  calculation  of  the  associated  Legendre  functions  (Gaussian 
quadratures,  alternate  series),  the  modulus  of  the  results  appeared  accurate 


15 


to  at  least  25  significant  figures.  Tabular  and  graphical  results  are 
presented  in  Appendix  A  and  Appendix  B,  respectively. 

In  routines  still  under  development  at  BRL,  special  series  and 

8Pm(cos0 ) 
v 

nonhomogeneous  recurrence  formulas  are  being  used  to  calculate  - jr - 

3Qv(cos0) 

and  - -  .  Also  being  programmed  are  asymptotic  expansions  to  extend 

the  range  of  |v|  to  very  large  values.  For  small  angles  (0<3O°),  Bessel 
function  expansions  will  be  used.8  Special  hypergeometric  series  along  with 
sines  and  cosines  will  be  used  for  larger  angles. 


8 

N.K .  Chu  Khrukidze,  "Asymptotic  Formulae  for  Legendre  Functions,  ”  USSR 
Computational  Mathematics  and  Mathematical  Physics ,  6,  1966 . 
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APPENDIX  A 

NUMERICAL  CALCULATIONS 
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CONICAL  LEGENDRE  FUNCTION  OF  THE  bECOhD  KIND 
COMPLEX  DEGREE 
OF  ORDER  M=  10 
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savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 


Name : 

Telephone  Number: 
Organization  Address: 


